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So far all ground-based retrievals of cloud properties assume that clouds vary only 
vertically, ignoring not only horizontal in-cloud structure but also broken cloudiness. 

This leads to misinterpretation of cloud properties and often makes their remote sensing 
impossible. We propose a new technique that retrieves cloud optical thickness for even 
broken clouds above green vegetation from surface measurements of zenith radiance in 
the red (RED) and near-IR (NIR) spectral regions. The idea of the method is simple and 
based on the very different spectral behavior of green vegetation and cloud liquid water 
drops. Since green vegetation reflects 30-50% of incoming radiation in the NIR and only 
5-10% in the RED region, ground measurements under thin clouds have little spectral 
contrast between RED and NIR, while thick clouds reflect much more of the surface- 
reflected radiation in the NIR than in RED. Based on this idea, we proposed to use points 
on the "RED versa NIR" plane for the retrieval of cloud optical thickness. 

The proposed retrieval method is applied to Cimel measurements at the Atmospheric 
Radiation Measurements (ARM) site in Oklahoma. Cimel, a multi-channel 
sunphotometer, is a part of AERONET - a ground-based network for monitoring aerosol 
optical properties. When conditions are inappropriate to make AERONET measurements 
that are suitable for aerosol studies, new measurements related to cloud physics are made 
instead. As such, several AERONET Cimel sunphotometers have been equipped with a 
new "cloud mode." This mode allows the Cimels to make measurements of zenith 
radiance when the Sun in blocked by clouds and thus to monitor cloud properties using 
the proposed retrieval method. 
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Abstract 

A new method for retrieving cloud optical depth from ground-based 
measurements of zenith radiance in the RED and near infrared (NIR) spectral regions is 
introduced. Because zenith radiance does not have a one-to-one relationship with optical 
depth, it is absolutely impossible to use a monochromatic retrieval. On the other side, 
algebraic combinations of spectral radiances such as NDCI (Marshak et al., 2000) while 
largely removing nouniquiness and the radiative effects of cloud inhomogeneity, can 
result in poor retrievals due to its insensitivity to cloud fraction. Instead, both RED and 
NIR radiances as points on the “RED vs. NIR” plane are proposed to be used for 
retrieval. The proposed retrieval method is applied to Cimel measurements at the 
Atmospheric Radiation Measurements (ARM) site in Oklahoma. Cimel, a multi-channel 
sunphotometer, is a part of AERONET - a ground-based network for monitoring aerosol 
optical properties. The results of retrieval are compared with the ones from Microwave 
Radiometer (MWR) and Multi-Filter Rotating Shadowband Radiometers (MFRSR) 
located next to Cimel at the ARM site. In addition, the performance of the retrieval 
method is assessed using a fractal model of cloud inhomogeneity and broken cloudiness. 
The preliminary results look very promising both theoretically and from measurements. 


1. Introduction 

The most common approach for retrieving cloud optical depth from ground-based 
observations uses downwelling fluxes measured by pyranometers in the 0.3- to 3.0- pm 
region of the solar spectrum (Leontieva and Stamnes, 1994; Boers, 1997). They are 
relatively cheap and included as standard equipment at many meteorological stations. In 
addition to broadband pyranometers, there are Multi-Filter Rotating Shadowband 
Radiometers (MFRSR) that infer the optical properties of clouds using downwelling 
fluxes measured at one or at several wavelengths in the visible and/or near infrared 
spectral regions (Min and Harrison, 1996; Leontieva and Stamnes, 1996). The key 
element in both retrieval techniques is the one-to-one mapping of the “observed” fluxes 
into cloud optical depth through (the use of) plane-parallel radiative transfer. Both 
methods are expected to work well only for completely overcast clouds (Ricchiazzi et al., 
1995) giving an effective optical depth for the whole sky since for inhomogeneous clouds 
each sky element contributes to the downwelling flux differently (Boers et al., 2000). To 
infer cloud optical depth locally one can assume to use a narrow-field-of view radiometer 



that measures radiances instead of fluxes (Pavloski and Ackerman, 1999). However, lack 
of one-to-one relationships between radiance and cloud optical depth (e.g., see Fig. 1 in 
Barker and Marshak, 2001 for zenith radiances), prevents the direct use of radiances 
measured either. 

Recently, Marshak et al. (2000) and Knyazikhin and Marshak (2000) proposed to 
exploit the sharp spectral contrast in vegetated surface reflectance across 0.7 pm 
wavelength to retrieve cloud properties from ground-based radiance measurements. The 
idea is to use ground zenith radiance measurements in two narrow spectral bands on each 
side of the sharp jump in vegetation albedo near 0.7 pm. Below 0.7 pm in the RED 
spectral region (648-680) nm), the chlorophyll in green vegetation absorbs 90-95% of 
solar radiation; thus its albedo is low. On contrast, above 0.7 pm in the near infrared 
(MR) spectral region, a green leaf reflects about 90% of incident radiation resulting in a 
high value of the albedo of vegetated surfaces. In the MR region, therefore, the green 
vegetation acts as a powerful reflector illuminating horizontally inhomogeneous clouds 
from below (see Fig. 5 in Marshak et al. (2000)). This provides the extra information 
needed to largely remove the ambiguity in measured downwelling radiance caused by 
radiative effects of the three-dimensional (3D) cloud structure. 

The spectral index proposed by Marshak et al. (2000) while substantially 
eliminates the effect of 3D cloud structure does not carry information on cloud fraction 
and, in general, can not correctly estimate the total amount of radiation reaching the 
surface. Thus in addition to zenith radiances, Barker and Marshak (2001) suggested to 
use downwelling spectral fluxes measured at the surface or from aircraft-based 
radiometers (Barker et al., 2002). Knowledge of downwelling fluxes substantially 
improves the retrieval though limits its applicability to sites equipped with MFRSR or 
other spectral flux measured radiometers (e.g. LI- 1800, ASD). 

Here we propose a new method for inferring cloud optical depth from spectral 
zenith measurements without requiring simultaneous measurements of spectral fluxes. 
Thus it can be used to monitor clouds from remotely located Cimel multi-channel 
sunphotometers that are the main part in of the AERONET - a ground-based monitoring 
network that consists of identical multi-channel radiometers for assessing aerosol optical 
properties and validating their satellite retrievals (Holben et al., 1998). 

The idea of the method comes from the “RED versa MR” method used to retrieve 
leaf area index (LAI) from satellite measurements (Knyazikhin. et al., 1998; Shabanov et 
al., 2002). In addition to cloud optical depth, the new retrieval method also infers a 
“radiatively effective” cloud fraction. Similar approach in remote sensing is a Nakajima- 
King technique for retrieving cloud optical depth and effective particle radius from 
(Nakajima and King, 1990). 

The plan of the paper is the following. We start from the physical background 
deduced from Cimel measurements under clear and cloudy conditions (section 2); then 
we discuss cloud index and its shortcomings (section 3). Section 4 introduces a new 
method and which will be applied to Cimel data (section 5). Surface reflective properties 



as a key element in the retrieval technique is discussed in section 6, while section 7 
validates the method using a fractal-based model of cloud inhomogeneity. Finally, 
section 8 summarizes the results. 


2. Three main atmospheric cases observed in ground zenith radiance measurements 

For proof-of-concept measurements (Wiscombe et al., 2000), Fig. 1 shows a 22 
min fragment of zenith radiance measured by a ground-based Cimel multi-channel 
sunphotometer pointed straight up. Cimel has a narrow field of view of 1.2° and four 
filters at 0.44, 0.67, 0.87 and 1.02 pm that are designed for retrieving aerosol properties 
in clear sky conditions. In our example, Cimel measured radiance at 20 sec. temporal 
resolution. 

There are three distinct regions in Fig. 1 : (from left to right) a single unbroken 
cloud, broken clouds, and a clear sky. For clear sky conditions, due to Rayleigh 
scattering and optically thicker aerosol at smaller wavelengths, zenith radiance increases 
with decrease in the wavelength from 1.02 to 0.44 pm. By contrast, for cloudy 
conditions radiances in channel 0.44 and 0.67 pm are almost indistinguishable; this also 
true for channels 0.87 and 1.02 pm. This is a clear indication that in the presence of 
clouds, the spectral contrast in surface albedo dominates over Rayleigh and aerosol 
effects. In contrast to the small fluctuations typical for clear and even cloudy skies, 
broken clouds show sharp changes in radiances around cloud edges. 

To be more formal, based on photon cloud-vegetation interactions we distinguish 
three main cases: 

(1) Atmosphere dominates. In this case, 

4.44 >:> 4.67 > 4.87 > 4.02 (l a ) 

and aerosol optical properties can be retrieved. 

(2) (Vegetated) surface and cloud dominates. In this case, 

•4.44 ~ 4.67 < 4.87 ~ 4.02 0^) 

and cloud optical properties can be retrieved (provided we know the surface albedo). 

(3) Transition between the first two case s, characterized by rapid changes between the 
“order" of 4. from cloudy to clear and back. In this case, neither aerosol nor cloud 
properties can be reliably retrieved using only one wavelength. However, as will be seen 
below, the two-wavelength retrieval for broken clouds can be almost as successful as for 
an overcast sky. 



3. Cloud Index and its shortcoming 


By analogy with Normalized Difference Vegetation Index (NDVI, Tucker, 1979; 
Verstraete and Pinty, 1996), Marshak et al. (2000) and Knyazikhin and Marshak (2000) 
proposed to use the Normalized Difference Cloud Index (NDCI) defined as a ratio 
between the difference and the sum of two normalized zenith radiances measured for two 
narrow spectral bands in the NIR (0.87 pm) and RED (0.67 pm) spectral regions, 

NDCI= Im zIm b. ( 2 ) 

T + J v ' 

1 NIR ^ 1 RED 

Compared to a two-valued optical depth versa zenith radiance relationship, that makes its 
retrieval absolutely impossible (see one-dimensional (ID) curves in Figs. 2a and b), the 
NDCI is a monotonic function with respect to optical depth (see ID curves in Fig. 2c). In 
contrast to any conventional method of estimating cloud optical depth from the surface 
that uses either broadband (Leontieva and Stamnes, 1994) or a single wavelength (Min 
and Harrison, 1996; Leontieva and Stamnes, 1996) and is expected to work well only for 
overcast clouds (Boers et al., 2000), the NDCI-based retrieval technique is much less 
sensitive to cloud structure. The sensitivity is weak because the NDCI-based method 
eliminates the part of downward radiation that did not have interactions with surface; this 
radiation is the most sensitive to both illumination conditions and cloud inhomogeneity 
(Marshak et al, 2000, Barker and Marshak, 2001). In addition, the NDCI is almost 
insensitive to the Solar Zenith Angle (SZA); consequently, optical depth of a cloud 
illuminated under SZA=80°, can be retrieved as good as the one illuminated under 
SZA=45°. This is a valuable feature of a retrieval method since most current techniques 
fail to perform reliable retrievals for large SZA. The NDCI-based approach first extracts 
radiation reflected by clouds and then performs retrieval; hence its weak sensitivity to the 
SZA (compare to Kaufmann et al., 2000 for NDVI). 

As follows from Eqs. (la)-(lb), the NDCI will be negative for a clear sky and 
positive for an overcast sky. In case of broken clouds, NDCI can take on either positive 
or negative values, depending whether there is a cloud in the zenith direction or not. 

The first shortcoming of the ID NDCI-based retrieval technique comes from the 
underestimation of zenith radiance for large optical depth in NIR. Indeed, Figs. 2a and 
2b show 3D and ID zenith radiances calculated for black (RED) and bright surface 
(NIR). We see that in the RED spectral region, 3D radiances are scattered around a 
theoretical ID curve, while in NIR, ID radiance systematically underestimates 3D 
radiances for large optical depths. This is understandable since for 3D clouds, more 
radiation is transmitted through; thus more radiation is reflected back from thick clouds 
to the surface. 

Another shortcoming of the NDCI is the same as for any spectral indices (e.g., 
Diner et al., 1999; Tian et al., 2000): the spectral information is reduced to one number 
by algebraic transformation. In other words, instead of two spectral values of zenith 



radiances in RED and NIR, we keep only one , NDCI. Indeed, each point on a RED versa 
NIR plane has two coordinates: 


V ~ + (nz* (^ a ) 

cc — arctan(/ gFn 1 1 m R ) . (3b) 


However, all points on a RED versa NIR plane with a similar a have the same 

NDCI = l ~ tzn J x (3c) 

1 + tana 

independent on rj. We will show below how this shortcoming can be overcome if both 
NIR and RED wavelength are used instead of the NDCI. 


4. RED versa NIR plane 


Any ground measurements of radiance, I, can be represented as a sum of two 
components (e.g., Box et al., 1988): the radiation calculated for a non-reflecting surface, 
7 0 , and the radiation due to surface-cloud interactions, 


/=/ o 


, PTJ, 

1 -pR 


(4) 


Here, R is the spherical albedo for isotropically illuminated clouds, 7^ is the transmittance 
for non-reflecting surface, and I 5 is the radiance generated by an isotropic source located 
at the surface with albedo p. 


Following Baker and Marshak (2001) we will define T 0 and R as 

To = T :lear ( 1 ~A C ) + T Aauiy A Q , (5a) 

R = RdenA 1 *^c) 4" Reload y^c, (^b) 

where A c is a cloud fraction. Assuming for simplicity that T daa = 1, T cioaiy = T^, /? dear = 0, 
and 7? cloudy == R^, where T’opp and /?pp are plane-parallel transmittance for non-reflecting 
surface and isotropic albedo, respectively. Thus, Eqs. (5a) -(5b) can be rewritten as 

F 0 = 1 -A c + ToppA,., (6a) 

R^R^Ao- (6b) 

Now we assume that in Eq. (4) only surface albedo p is wavelength dependent while 7 0 , 

7 S , T 0 and R are functions of only cloud optical depth x. It follows from these assumptions 
and from (6a) -(6b) that 



(7a) 


and 


= W+ 1 


1 Pi2£D^C^»p(^) 




1-P«a4*-(T) 


(7b) 


Let us consider Eq. (6a). It is easy to show that 


if 


O^T^T^l ( 8 ) 

0<A C <1. (9) 


The middle inequality in (8) expresses a common knowledge that transmittance for 
inhomogeneous (or broken) clouds is larger than its plane-parallel counterpart. However, 
it is also known that for broken clouds transmittance can exceed unity because of photons 
scattered downwards from cloud edges. Thus the last inequality in (8) is not necessarily 
hold for broken clouds. To illustrate this, Fig. 3 shows broadband measurements of 
downward flux for two days in August 1997. We can clearly see that downward flux for 
broken clouds (between 9 and 10 a.m. on Aug. 17 and around noon on Aug. 15) exceeds 
the one for clear sky thus producing transmittance larger than 1. More accurately, at 
12:05 and SZA=19.2 on Aug. 15 the transmittance even exceeds 1.2. Can this effect be 
simulated in the frame of plane-parallel radiative transfer? 


Indeed, if one releases condition (9) and allows “cloud fraction” to be negative 
then, as it follows from Eq. (6a), T 0 will exceed 1 and be able to simulate reflectance 
from cloud edges. Not surprisingly that in this case, T 0 as a function of optical depth 
looks more like reflectance than transmittance (convex vs. concave). Figure 4 shows T 0 
as a function of x and A c . Note that A e = 0 gives T 0 = 1 while negative A c correspond to T 0 
> 1 . 


To illustrate this behavior Fig. 5a shows an 8-km fragment of a simulated broken 
cloud field. As expected, the cloud edge around 38.4 km looks much brighter than the 
rest of the cloud for both RED and NIR radiation. However, when one moves away from 
the cloud edge, / ffiD decays much faster than / MR . Only about 300 m away from the edge, 
/red has been already stabilized while / NIR is still strongly affected by the illumination by 
the bright surface. (Note from the ID curves in Figs. 2a and 2b that max(/ RED ) = 0.53 
while max^REp) = 0.62). No positive ^4 c is able to model this effect in the frame of plane- 
parallel theory. Indeed, panel 5b shows a 2D Look-Up-Table (LUT) of / NIR versa /red as 
functions of x and A c described by Eqs. (7a) - (7b). We see that the data-points that are 
above the A e = 0 curve need negative ^ c to be matched by ID calculations! Obviously, 
these data-points correspond to transmittance T 0 > 1 . Also note that all cloudy data- 
points as well as the LUT itself is above the diagonal I mR = /re D , i.e. 7^ > /re D . The same 
data-points and LUT will occupy the whole plane if the difference 4® - /re D is plotted 
against /re D or the sum / NIR + / RED . The later one (panel 5c) is preferable since any line 
through die origin point (0,0) corresponds to NDCI defined in Eq. (2). 



As an example, Fig. 6a shows a NIR-RED versa NIR+RED plane with real Cimel 
measurements at the Atmospheric Radiation Measurements (ARM) site in Oklahoma 
(97.48° W, 36.61° N) on My 28, 2002. For SZA = 62° and surface albedos Pre D =0.092 
and p N]R =0.289, it illustrates a DISORT calculated LUT with and / RED as functions of 
x and A c , and three groups of data-points (10 points per a group, see next section for 
explanations) measured by Cimel. The data-points, while having different values of 
and / RED thus being located at different positions on the plane, have almost the same 
NDCI (a straight line through the point (0,0)); hence the same optical depth x (80 for A c = 
1 .0). However, as it follows from the LUT, these groups correspond to different pairs: 
(A = 0.9; x=28), (4 C =0.8; x=22), and (4 C =0.4; x=12); thus have different optical depths. 
Figure 6b shows a Total Sky Image (TSI) taken at 14:00 UMT close to the middle cluster 
of points in panel 6a. Note that A c = 0.8 is not a “real” cloud fraction but rather a 
“radiatively effective” one in the sense that ID radiative transfer gives the same values of 
Ink and I mD as the measured (3D) ones. 

The LUTs shown in Fig. 6 and Figs. 5b and 5c for retrieval cloud optical depth 
and cloud fraction can be viewed as a ground-based counterpart of the Nakajima-King 
retrieval of cloud optical depth and effective particle radius from reflected solar radiation 
(Nakajima and King, 1990). In the Nakajima-King retrieval, visible channel is more 
sensitive to cloud optical depth and less sensitive to a particle size while the absorbing 
NIR channel is more sensitive to a particle size than to optical depth. Similar, for the 
ground-based retrieval above vegetation, the RED channel is much more sensitive to 
optical depth while the NIR is more sensitive to cloud fraction. Next section will discuss 
the retrievals in more details. 


5. Retrieval from Cimel measurements 

Cimel spectrometers are the main instruments in the AERONET - a ground based 
monitoring network that consists of identical multi-channel radiometers for assessing 
aerosol optical properties and validating their satellite retrievals (Holben et al., 1998). 
Based on the above ideas, we have developed a method for Cimels to monitor cloud 
optical depth by using “idle time” inappropriate for aerosol studies for taking 
measurements of zenith radiance at 0.44, 0.67, 0.87 and 1.02 pm. Every 13 min, if the 
Sim is blocked by clouds, the Cimel points straight up and takes 10 measurements with a 
9 sec. time interval. 

The main idea of our method is to retrieve cloud optical depth using both 0.67 pm 
(RED) and 0.87 pm (NIR) zenith radiances; surface albedos p^ and p NYR are assumed to 
be known from independent measurements. Other two channels (0.44 and 1.02 pm) are 
used to select cases where the spectral contrast in surface albedo dominates over 
Rayleigh and aerosol effects (i.e., when Eq. (lb) is hold) and the proposed method will 
likely work. Another consistency check is related to Cimel’ s second collimator that has 
the same field of view of 1.2° but 10 times larger aperture; thus it sees different cloud 
pieces. 



Figure 7a illustrates the retrieval of cloud optical depth from Cimel measurements 
using the above method. The retrieved optical depths are also compared with the ones 
from MicroWave Radiometer (MWR) (assuming a constant droplet effective radius of 1 0 
pm) and from Multi-Filter Rotating Shadowband Radiometer (MFRSR) (courtesy of Q. 

L. Min; for details see Min and Harrison, 1996). For this overcast day, generally we see 
very good agreement, though there are a few discrepancies. For the most homogeneous 
plane-parallel cloud around 14:30 UMT (panel 7b 1) all three methods show similar 
optical depth x around 20; the LUT in Fig. 7b2 estimates the radiatively effective cloud 
fraction A c as 0.8. For more inhomogeneous though still overcast clouds as the one at 
17:50 (panel 7cl) the agreement is not as good as for the homogeneous cloud since 
MFRSR retrieves a more effective optical depth (in its field of view) than our algorithm 
does. Figure 7c2 shows the LUT with scattered Cimel dots; the most inhomogeneous 
part (at 17:52) corresponds to the largest difference between and 1^ that estimates A c 

~ 0.5 and x ~ 90. 

Another example (Fig. 8a) corresponds to broken clouds (Figs. 8b and 8c). There 
are two remarkable features highlighted in this day. First, the optical depth is retrieved 
for SZA=81° (at 12:30-12:40, see panel (b)); this is not very common for any (satellite or 
ground-based) optical depth retrievals. Second, it is retrieved for very broken-cloud 
conditions; the cloud optical depth for 30 min. jumps from almost 0 (17:52) to 90 (1 8:08) 
and then backs to clear sky (18:22). This is also shown in Fig. 8d where we see the 
rapidly changing order of Cimel zenith radiances as in Fig. 1 and Eqs. (la) - (lb). As 
expected, the optical depth field retrieved from MFRSR around 18:00 is much smoother 
and represents an “effective” optical depth for the sky shown in panel (c). 


6 . Surface 

Since the retrieval algorithm is based on surface-cloud interaction, knowledge of 
surface reflectance and its heterogeneity around Cimel is absolutely crucial for the 
success of the algorithm. 

To estimate surface albedo we use Moderate Resolution Imaging 
Spectroradiometer (MODIS) data averaged over 16 days (Schaaf et al., 2002). Figure 9 
illustrates surface albedo for 0.648 and 0.858 |im around the ARM site for July 28, 2002. 
The average values are about 0.09 for 0.648 (im and 0.29 for 0.858 pm; they are 
consistent with Multi-angle Imaging Spectro-Radiometer (MISR) data (REFERENCE) 
retrieved on July 24, 2002 (p 0 . 672 = 0-12, p 0 . 867 = 0.28; averaged over (27.5 km) 2 : <p 0 . 672 > = 
0.10 with standard deviation (sdev) = 0.02 and <p 0 . 867 > = 0-28 with sdev=0.03. Pavloski 
et al. (2002) used MISR surface albedo data to retrieve cloud optical depth from 
simultaneous measurements of downward radiance and flux at the ARM site in 
Oklahoma applying algorithm proposed by Barker and Marshak (2001). 


For the theoretical study based on cloud fields inferred from Landsat imagery, 
Baker and Marshak, (2001) found that as long as the uncertainties in surface albedo have 
the same sign (both are either overestimated or underestimated), the algorithm performs 



well. Furthermore, when the NIR albedo is overestimated but the RED albedo is 
underestimated, errors in the retrieved optical depth are not severe. However, in the 
opposite case, the algorithm underestimates multiple reflectance in the bright band and 
greatly overestimates optical depths. A detailed analysis of the sensitivity to surface 
albedo can be found in Beaulne et al. (2003). Similar tendencies are also valid for the 
present algorithm. 

Surface reflectance at the ARM site in Oklahoma is very heterogeneous (Li et al., 
2002). Figure 10 is based on Solar Spectral Flux Radiometer (SSFR) data (Pilewskie et 
al., 2003) reported in Michalsky et al. (2003). It shows surface albedo in RED and NIR 
spectral regions measured on flight legs covering the area of about 25 km around the 
ARM site (panel (a)) on April 5, 2000. In addition, data for the nearest 0.0001 radian 
(about 600 m) to the site was also selected. We see that the range of surface albedos in 
the area of 1 km around Cimel on that day was from 0.056 to 0.1 18 for the RED and from 
0.325 to 0.434 for the NIR with = 0.085 ±0.013 and = 0.375 ±0.021. These 
average results are not necessarily match surface albedo measured from 10-m and 25-m 
towers located exactly at the ARM site over wheat and (dead) grass (Li et al., 2002, 
Michalsky et al., 2003). However, the big range in surface albedos does not prevent the 
use of the proposed method since the contribution from surface is integrated over several 
km (depending on cloud base height) and the local surface properties are less important 
than its average or “effective” characteristics (Knyazikhin and Marshak, 2000). It can be 
shown that in spite of strong surface heterogeneity the value of surface reflectivity around 
the ARM site stabilizes after averaging over 1 - 2 km for both RED and NIR spectral 
regions (Wiscombe et al., 2000). 


7. Validation 

In addition to validation with retrievals from other ground-based measurements 
(MWR and MFRSR) that are illustrated in Figs. 7 and 8, we show a model-based 
validation where the true value of optical depth is known. We used 10 realizations of a 
10-cascade fractal cloud model of broken clouds from Cahalan et al., 1994 and Marshak 
et al., 1998 (see Fig. 2 captions for details). The cloud fraction was 81% and the total 
number of nonzero cloud optical thicknesses was (0.81 x 2 10 x 10 = 8300). The average 
in-cloud T was 13 with standard deviation of 10.5. For simplicity, we had “black” and 
“bright” surface with albedos = 0.0 and p^ = 0.5, respectively. The domain 
averaged retrieved values for oblique illumination of SZA = 60° and conservative 
scattering gave very close results of mean 13.1 and standard deviation 10.8. 

Figure 11a shows a scatter plot of retrieved versa true optical depth on a pixel-by- 
pixel basis for a pixel size of 25 m. We see that while the “NDCF’ method overestimates 
cloud optical depth values especially for t > 20 and shows a lot of saturated values of t = 
1 00, the “RED vs. NIR” method gives on average unbiased estimates. When averaged 
over 8 pixels (200-m resolution pixels), the retrieval substantially improves (Fig. 1 lb); 
the average absolute error in x is 1.8 down from 2.4 (3.0 for in-cloud pixels) for 25-m 
resolution pixels. 



Finally, Figs. 12a and b illustrate histograms of cloud optical depth (true and 
retrieved) on a lin-lin and log-line scale, respectively. We can clearly see that the “RED 
vs. NIR” method correctly retrieves the distribution of X for both small and large values. 
Analyzing the cumulative histogram of absolute errors (Fig. 13) we can conclude that in 
50% cases the absolute error in the retrieval of optical depth is smaller than 1, in 75% 
cases smaller than 3, and in 90% cases is smaller than 5. In contrast, the “NDCI” method 
strongly overestimates the retrieved values for at least x > 40; this makes it completely 
unacceptable for the retrieval of large x. 


8. Summary 

Ground-based remote sensing of cloud optical properties with broadband 
pyranometers (e.g., Leontieva and Stamnes, 1994) and/or narrowband radiometers (Min 
and Harrison, 1996) that measure downwelling fluxes is limited to overcast conditions 
(Boers et al., 2000). Retrieved values of cloud optical depth have an effective rather a 
local character (Barker and Marshak, 2001). Monochromatic narrow field-of-view 
measurements are not widely applicable for retrieval because cloud optical depth is a 
two-value function of zenith radiance. Among others these shortcomings prevent us from 
creating a ground-based cloud properties monitoring network as the AERONET for 
assessing aerosol properties (Holben et al., 1998). 

To remove the ambiguity in measured downwelling radiance, it was suggested to 
exploit a spectral contrast in surface properties between RED and NIR spectral regions 
(Marshak et al., 2000). An algebraic combination of NIR and RED zenith measurements 
or, the so-called Normalized Difference Cloud Index (NDCI), largely removes not only 
the double valued relationship between radiance and optical depth but also the radiative 
effects of the 3D cloud structure. However for accurate retrievals it requires additional 
information on downwelling fluxes (Barker and Marshak, 2001 and Barker et al., 2002) 
and the NDCI (or other similar indices) along cannot be used for broken cloud 
conditions. 

Here we proposed a new method that overcomes the limitations of the NDCI; 
instead of using one value of an algebraic combination between two wavelengths, we use 
both wavelengths independently in the NIR vs. RED plane. Mapping NIR and RED 
measurements of zenith radiance into a two-dimensional plane-parallel radiative transfer 
based look-up-table (LUT) allowed us not only to retrieve cloud optical depth x but also 
to estimate cloud fraction A c . However, since plane-parallel radiative transfer is unable to 
fit all 3D measurements (for example, leakage from cloud edges results in cloud 
transmittance larger than unity), instead of a “real” cloud fraction (the ratio of cloudy 
pixels to the total number of pixels) we retrieve a radiatively effective cloud fraction (or 
index of local cloud inhomogeneity) that matches 3D measurements with ID 
calculations. This does not preclude it from to be negative. The retrieved values A Q are 
not used per se but as a degree of cloud inhomogeneity for the retrieval of cloud optical 
depth. 



This method assumes that surface albedo in both RED and NIR spectral regions is 
known. Because of the proposed global character of the cloud monitoring network, the 
high resolution satellite data is the only available source of surface properties. As an 
example, here we applied MODIS surface albedo data averaged over 16 days. Even in 
case of heterogeneous surface as the one around the ARM site in Oklahoma (Li et al., 
2002; Michalsky et al., 2002), we found that the proposed method can be successfully 
used since the local surface properties are less important than their average characteristics 
(Knyazikhin and Marshak, 2000). 

In addition to the comparison with MFRSR (courtesy of Qilong Min) and MWR 
retrievals (the later needs an assumption on the drop size) that are located next to Cimel 
at the ARM site, we used realizations of a fractal model of cloud inhomogeneity (Cahalan 
et al., 1994) and of broken clouds (Marshak et al., 1998). Since the “true’ optical depth 
was known, we were able to validate the performance of the proposed method. 75 % of 
the retrieved optical depth out of about 10000 pixels with 25 m resolution had absolute 
errors smaller than 3. Averaging over 200 m substantially improves the retrieval. 

To conclude, the preliminary results look very promising both theoretically and 
from measurements. As soon as the method matures, it will be systematically applied to 
AERONET “cloud” measurements that are inappropriate for aerosol studies to monitor 
cloud optical depth. 
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Figure Captions 

Figure 1. Zenith radiance measured by a Cimel sunphotometer at Greenbelt, MD on May 
24, 1999. Four channels 0.44, 0.67, 0.87 and 1.02 pm are used. The measured radiance is 
normalized by the solar flux at the TOA in the corresponding spectral interval. 

Figure 2. Downward radiances calculated by Monte Carlo methods for “black” surface 
(RED, p = 0.0, panel (a)), and “bright” surface (NIR, p = 0.5, panel (b)). Pixel size is 25 
m, SZA = 60°, G5 0 = 1.0, Henyey-Greenstein scattering phase function. Horizontal 
distribution of cloud optical depth is simulated by a 10-steps bounded cascade model 
(Cahalan et al., 1994) with parameters <x> - 13, P = 1.4 and p = 0.35. The average 
geometrical cloud thickness is 300 m; cloud base height is 1 km. Holes are added as in 
Marshak et al. (1998). The results of ID radiative transfer calculations from DISORT 
(Stamnes et al., 1988) are added for convenience, (c) The same as in panels (a) and (b) 
but for the NDCI defined by Eq. (2). 

Figure 3. Ground based measurements of downward radiation at La Jolla, CA in August, 
1997. Left panel is typical for La Jolla: early morning overcast sky, then broken 
cloudiness that clears up by noon. Right panel represents an almost entirely overcast day 
with a few scattered clouds around noon. Solar zenith angles are added for clarity. 

Figure 4. Transmittance T 0 as a function of optical depth x and cloud fraction A c . 
Horizontal axis consists of 15 ranges of x from 0 to 100. Each x-range corresponds to 
different A c from A c = 1 to A c = -0.4 including the case of^ c = 0 which corresponds to T 0 
= 1 . 

Figure 5, Zenith radiances 7,® and 1 ^ for the same model as in Fig. 2. (a) An 8-km 
fragment of / RED , / NIR (top panel, right vertical axis), and optical depth x (bottom panel, 
left vertical axis) versa horizontal position x (km), (b) Zenith radiances for RED and NIR 
as functions of optical depth x and cloud fraction A c [Eqs. (7a) and (7b)] calculated using 
DISORT (Stamnes et al., 1988); x changes from 0.5 to 100 while A c changes from 0.0 to 
1 .0. Note that both and 7^ are double-value functions with respect to x. Dots 
correspond to data-points from panel (a) between 38.8 (I MR < 0.8) and 44 km (208 dots), 
(c) The same as in panel (b) but for NIR-RED versa NIR+RED. 

Figure 6. NIR-RED versa NIR+RED plane with Cimel radiances measured at the ARM 
Oklahoma site on July 28, 2002. (a) SZA = 62° ± 3° and surface albedos p RED =0.092 and 
p N1R =0.289. /red and / NIR are calculated using DISORT. Optical depth x changes from 0.5 
to 80 while cloud fraction A c changes from 0.2 to 1.0. 30 dots (10 for each cluster) 
correspond to Cimel measurements taken around 13:45, 13:58, and 14: 1 1 UMT, 
respectively (towards decreasing NIR+RED). A straight line through (0,0) is the NDCI = 
0.08. (b) a Total Sky Image taken at 14:00 UMT with SZA=62.75°. 



Figure 7. Retrievals for July 3, 2002 at the ARM Oklahoma site, (a) retrievals based on 
MWR (assuming droplet effective radius of 10 pm), Cimel radiances (using reflectances 
from a vegetated surface with albedos Pred = 0- 119 and p Nm =0.302), and MFRSR fluxes 
(using algorithm proposed in Min and Harrison, 1996). (bl) a Total Sky Image for SZA 
= 52.3° taken at 14:36 UMT. (b2) LUT for SZA = 53° ± 1°. Optical depth x changes from 
0.5 to 1 00 with increment of 2.5 while cloud fraction A c changes from -0.6 to 1 .0 with 
increment of 0.2. 20 dots correspond to Cimel measurements taken between 14:30-14:37 
UMT. (cl) the same as in (bl) but for SZA = 16.3° taken at 14:50 UMT. (c2) the same as 
in (b2) but for SZA = 16° ± 1°. 20 Cimel measurements were taken between 17:43-17:52 
UMT. 

Figure 8. Retrievals for August 8, 2002 at the ARM Oklahoma site, (a) retrievals based 
on MWR (assuming droplet effective radius of 10 pm) and Cimel radiances (surface 
albedos are Pre D = 0-096 and p m =0.338). (b) a Total Sky Image for SZA = 79.6° taken at 
12:40 UMT. (c) the same as in (b) but for SZA = 21.9° taken at 18:00 UMT. (d) 30 min. 
of normalized Cimel zenith radiances at four channels: 0.44, 0.67, 0.87 and 1.02 pm. 

Figure 9. Surface albedo retrieved from MODIS on July 28, 2002. 25 by 25 km 2 area 
around the ARM site in Oklahoma; its location is at (0,0). (a) 0.648 pm, at the ARM site 
Po.648 = 0.092, averaged over (25 km) 2 <p 0 .648 > = 0.091 with standard deviation (sdev) = 
0.016. (b) 0.848 pm; at the ARM site p 0858 = 0.289, averaged over (25 km) 2 <p 0 . 8 5 8 > = 
0.292 with sdev = 0.043. 

Figure 10. Surface around the ARM site in Oklahoma, (a) Area around the ARM site 
with a daisy pattern flown by the Twin Otter aircraft on April 5, 2000 (courtesy of 
Alexander Trishchenko). (b) Surface albedo at 0.86 pm versa 0.66 pm measured from the 
aircraft by SSFR on April 5, 2000 (for details see Michalsky et al., 2003). Small gray dots 
correspond to all data collected around the area shown in panel (a) (roughly 10-15 km 
from the ARM site). Big black dots are the closest to the ARM site measurements 
(several middle points) taken within 0.0001 radian from the site (around 600 m). For 
them, <p 0 .66> = 0.085, sdev = 0.013 and <p 086 > = 0.375, sdev = 0.021. 

Figure 11. Scatter plot of retrieved versa true cloud optical depth for the same bounded 
cascade model as in Fig. 2 (Cahalan et al., 1994, Marshak et al., 1998). Cloud fraction is 
81%. Ten realizations of bounded cascades are used. SZA = 60°, Two retrieval methods 
are compared: the “NDCI” and the “RED vs. NIR”. (a) original resolution of 25 m. (b) 
averaged over 200 m. 

Figure 12. Histogram of cloud optical depth x from Fig. 11a. Pixel resolution is 25 m. (a) 
lin-lin plot, (b) log-lin plot. 

Figure 13. Cumulative histogram of absolute errors in the retrieval of cloud optical depth 
from Fig. 1 la. Pixel resolution is 25 m. 
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Figure 1 . Zenith radiance measured by a Cimel sunphotometer at Greenbelt, MD on 
May 24, 1999. Four channels 0.44, 0.67, 0.87 and 1.02 pm are used. The measured 
radiance is normalized by the solar flux at the TOA in the corresponding spectral interval. 
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Figure 2. Downward radiances calculated by Monte Carlo methods for “black” surface 
(RED, p = 0.0, panel (a)), and “bright” surface (NIR, p = 0.5, panel (b)). Pixel size is 25 
m, SZA = 60°, GJ 0 = 1.0, Henyey-Greenstein scattering phase function. Horizontal 
distribution of cloud optical depth is simulated by a 10-steps bounded cascade model 
(Cahalan et al., 1994) with parameters <T> = 13, P = 1.4 and p = 0.35. The average 
geometrical cloud thickness is 300 m; cloud base height is 1 km. Gaps are added as in 
Marshak et al. (1998). The results of ID radiative transfer calculations from DISORT 
(Stamnes et al., 1988) are added for convenience, (c) The same as in panels (a) and (b) 
but for the NDCI defined by Eq. (2). 
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Figure 3. Ground based measurements of downward radiation at La Jolla, CA in August, 
1 997. Left panel is typical for La Jolla: early morning overcast sky, then broken 
cloudiness that clears up by noon. Right panel represents an almost entirely overcast day 
with a few scattered clouds around noon. Solar zenith angles are added for clarity. 
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Figure 4. Transmittance T 0 as a function of optical depth x and cloud fraction A c . 
Horizontal axis consists of 15 ranges of x from 0 to 100. Each x-range corresponds to 
different A c from A c = 1 to A c = -0.4 including the case of A c = 0 which corresponds to T 0 
= 1 . 
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Figure 5. Zenith radiances 1^ and 7^ for the same model as in Fig. 2. (a) An 8-km 
fragment of 1^, (top panel, right vertical axis), and optical depth x (bottom panel, 
left vertical axis) versa horizontal position x (km), (b) Zenith radiances for RED and NIR 
as functions of optical depth x and cloud fraction A c [Eqs. (7a) and (7b)] calculated using 
DISORT (Stamnes et al., 1988); x changes from 0.5 to 100 while A c changes from 0.0 to 
1.0. Note that both and 7 Nm are double-valued functions with respect to x. Dots 
correspond to data-points from panel (a) between 38.8 (7,^ < 0.8) and 44 km (208 dots), 
(c) The same as in panel (b) but for 7 MR - 7^ versa 7™ + 7™. 










Figure 6. NIR-RED versa NIR+RED plane with Cimel radiances measured at the ARM 
Oklahoma site on July 28, 2002. (a) SZA = 62 0 ± 3° and surface albedos p RED =0.092 and 
p NIR =0.289. /red and / NrR are calculated using DISORT. Optical depth x changes from 0.5 
to 80 while cloud fraction A c changes from 0.2 to 1.0. 30 dots (10 for each cluster) 
correspond to Cimel measurements taken around 13:45, 13:58, and 14:1 1 UMT, 
respectively (towards decreasing NIR+RED). A straight line through (0,0) is the NDCI « 
0.08. (b) a Total Sky Image taken at 14:00 UMT with SZA=62.75°. 
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Figure 7. Retrievals for July 3, 2002 at the ARM Oklahoma site, (a) retrievals based on 
MWR (assuming droplet effective radius of 10 pm), Cimel radiances (using reflectances 
from a vegetated surface with albedos p^O.l 19 and p Nm =0.302), and MFRSR fluxes 
(using algorithm proposed in Min and Harrison, 1996). Solar zenith angle is added for 
convenience, (bl) a Total Sky Image for SZA = 52.3° taken at 14:36 UMT. (b2) LUT for 
SZA = 53° ± 1°. Optical depth x changes from 0.5 to 100 with increment of 2.5 while 
cloud fraction A c changes from -0.6 to 1.0 with increment of 0.2. 20 dots correspond to 
Cimel measurements taken between 14:30-14:37 UMT. (cl) the same as in (bl) but for 
SZA = 16.3° taken at 17:50 UMT. (c2) the same as in (b2) but for SZA = 16° ± 1°. 20 
Cimel measurements were taken between 17:43-17:52 UMT. 
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Figure 8. Retrievals for August 8, 2002 at the ARM Oklahoma site, (a) retrievals based 
on MWR (assuming droplet effective radius of 10 pm), Cimel radiances (surface albedos 
are Pre^O. 096 and p NIR =0.338), and MFRSR (Min and Harrison, 1996). Solar zenith 
angle is added for convenience, (b) a Total Sky Image for SZA = 79.6° taken at 12:40 
UMT. (c) the same as in (b) but for SZA = 21.9° taken at 18:00 UMT. (d) 30 min. of 
normalized Cimel zenith radiances at four channels: 0.44, 0.67, 0.87 and 1.02 pm. 
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Figure 9. Surface albedo retrieved from MODIS on July 28, 2002. 25 by 25 km 2 area 
around the ARM site in Oklahoma; its location is at (0,0). (a) 0.648 pm, at the ARM site 
Po. 648 = 0.092, averaged over (25 km) 2 <p 0 . 648 > = 0.091 with standard deviation (sdev) = 
0.016. (b) 0.848 pm; at the ARM site p 08J8 = 0.289, averaged over (25 km) 2 <p 0 . 858 > = 
0.292 with sdev = 0.043. 

(for MISR on July 24, 2002: p 0 . 672 = 0.123, averaged over (27.5 km) 2 <p 0 . 672 > = 0.103 with 
sdev = 0.017. (b) 0.867 pm; at the ARM site p 0 . 867 = 0.278, averaged over (27.5 km) 2 
< Po. 867 > = 0.276 with sdev = 0.030) 
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Figure 10. Surface around the ARM site in Oklahoma, (a) Area around the ARM site 
with a daisy pattern flown by the Twin Otter aircraft on April 5, 2000 (courtesy of 
Alexander Trishchenko). (b) Surface albedo at 0.86 |0.m versa 0.66 pm measured from 
the aircraft by SSFR on April 5, 2000 (for details see Michalsky et al., 2003). Small gray 
dots correspond to all data collected around the area shown in panel (a) (roughly 10-15 
km from the ARM site). Big black dots are the closest to the ARM site measurements 
(several middle points) taken within 0.0001 radian from the site (around 600 m). For 
them, <p 0 .66 > = 0 . 085 , sdev = 0.013 and <p 086 > = 0.375, sdev = 0.021 . 
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Figure 11 . Scatter plot of retrieved versa true cloud optical depth for the same bounded 
cascade model as in Fig. 2 (Cahalan et al., 1994, Marshak et al., 1998). Cloud fraction is 
81%. Ten realizations of bounded cascades are used. SZA = 60°. Two retrieval methods 
are compared: the << NDCF’ and the “RED vs. NIR”. (a) original resolution of 25 m. (b) 
averaged over 200 m. 
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Figure 12. Histogram of cloud optical depth x from Fig. 12a. Pixel resolution is 25 m. 
(a) lin-lin plot, (b) log-lin plot, (c) 
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Figure 13. Cumulative histogram of absolute errors in the retrieval of cloud optical depth 
from Fig. 11a. Pixel resolution is 25 m. 


